{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 2, demo 3\n",
    "\n",
    "Bayesian Data Analysis, 3rd ed\n",
    "\n",
    "\n",
    "Probability of a girl birth given placenta previa (BDA3 p. 37).\n",
    "Simulate samples from Beta(438,544), draw a histogram with quantiles, and do the same for a transformed variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import preliz as pz\n",
    "\n",
    "pz.style.use(\"preliz-doc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample from posterior Beta(438,544).\n",
    "# Obtain all draws at once and store them in vector 'theta'\n",
    "n = 10000\n",
    "theta = pz.Beta(438, 544).rvs(size=n)\n",
    "\n",
    "# Compute odds ratio for all draws\n",
    "phi = (1-theta)/theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/MAAAGbCAYAAACIxMC9AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOJ9JREFUeJzt3Xl4FGW6/vG7syckcABBkAEUoYOIKCiLyA9FFDyAuK8oOIiKYlzmBESdQR3x6IyoyOKIiIA6KCDMwAgOiAdZlUSWUUSCLLLEhU0hmD1dvz+YtGQS7HSn0/W+yfdzXVxX2Z3ufuq17656ut6q9jiO4wgAAAAAAFgjyu0CAAAAAABAcGjmAQAAAACwDM08AAAAAACWoZkHAAAAAMAyNPMAAAAAAFiGZh4AAAAAAMvQzAMAAAAAYBmaeQAAAAAALEMzDwAAAACAZWjmAQAAAACwDM08AAAAAACWoZkHapkffvhBTz/9tC699FJ16NBBffr00WuvvSafz+d2aQAAIAgjRoxQamqqli9f7nYpAFxAMw/UIuvWrdOAAQP09ttvy3EctWzZUrt379YLL7ygZ555xu3yAABAELKysiRJXq/X5UoAuIFmHqgltm3bpuHDh6uwsFAvvfSSli9frn/84x967bXXFBUVpVmzZmnXrl1ulwkAACohNzdX+/btU3Jyspo1a+Z2OQBcQDMP1AI+n0+PPPKIcnNz9dRTT6lfv37++y6++GL17NlTPp9Py5Ytc7FKAABQWdu2bZPjOByVB2oxmnmgFli8eLG2bNmic845R1dffXW5+8877zxJ0o4dOyJbGAAACMm2bdskSampqS5XAsAtMW4XAKD6zZgxQ5J02223VXh/cnKyJOmnn36KUEUAACBYBw4c0GuvvaaPPvpI33//vSRpzpw5Wr16ta699lrdfffdiolh9x6oLTgyD9Rwe/bs0RdffKG4uDhddtllFf5NUVGRJCkuLi6SpQEAgEr67LPPdOWVV+rNN9+Uz+fzN+3NmzfX3r179fLLL+sPf/iDy1UCiCSaeaCGW7lypSTpggsu8B+B/09Hjx6VJDVs2DBidQEAgMrZsWOHhg8friNHjmjMmDH66KOPFBUVpejoaC1cuFAvvfSSJGn+/Pn+6fcAaj7m4QA13IYNGyRJa9euDXheXcuWLSNREgAACMKTTz6pnJwcDR06VIMGDdLOnTuVl5en1q1bKz4+Xv369dO0adO0efNmrVmzhoviAbUEzTxQw3311VeSpLZt2yopKanCv/n8889VXFystm3bRrI0AAAQwL/+9S9lZGTolFNO0QMPPCBJ2rp1qySV2W63a9dOmzdv1uHDh12pE0Dk0cwDNZjjONq3b58kadKkSWrevHm5vzl8+LC6d++umJgYnXvuuZEuEQAA/Ip//vOfkqRLLrlEiYmJkn75ov6ss87y/53P55MkxcfHR7hCAG7hnHmgBjt27JgKCwslSY0aNarwb1avXi3HcdSlSxf/TgIAADDD+vXrJUk9evTw31Z6ZP7EZn737t2SpBYtWkSwOgBuopkHarDSq9R7PB4lJCRU+DeLFy+WJPXv3z9idQEAgMopbdJPP/10/20nnkInSTk5Odq0aZMkqXPnzhGtD4B7aOaBGqxu3bqKioqS4zjKy8srd/+uXbu0YsUKNWrUSFdeeaULFQIAgF9Tuv32eDySjp8ed+DAATVq1Mj/KzRLlixRUVGROnbsqKZNm7pWK4DIopkHarCYmBj/Feo///zzMvf5fD798Y9/lM/n04MPPsg5dgAAGKi0YS+dWr9lyxZJv0yxz8nJ0fjx4yVJd911V+QLBOAamnmghuvbt68k6cUXX9SxY8ckHT+XftSoUVq7dq169eqlG264wc0SAQDASXTv3l2SNHXqVB05cqTM+fKHDh1SWlqaDhw4oAEDBqh3795ulgogwjyO4zhuFwGg+hw5ckQDBgzQ/v37lZKSohYtWmjXrl3Kzc1Vjx49NHny5JOeTw8AANyVnZ2ta6+9Vj/99JNSUlIUHx+vgwcPqmnTpjp8+LAKCgrUp08fvfDCC4qLi3O7XAARRDMP1AL79u3T888/r08//VSFhYVq3bq1brzxRl1//fX+c/AAAICZduzYoeeff16ZmZn+WXZJSUlq27atbrnlFl155ZVsz4FaiGYeAAAAsEBBQYHOO+88JSQkaMOGDTTwQC3HOfMAAACABXbs2CGfz6czzzyTRh4AzTwAAABgg6+//lqS1KZNG5crAWACmnkAAADAAtu3b5cktW7d2uVKAJiAZh4AAACwQOmReZp5ABIXwAMAAAAAwDocmQcAAAAAwDKVauYdx1Fubq44iA8AgH3YjgMAUPNUqpnPy8vTkCFDlJeXV931wEBOfr6O3n6Ljt5+i5z8fLfLARAkMgy24+Ygj4BdyCxMxjR7AAAAAAAsE+N2AbBAVJSi257lXwZgGTIMmIM8AnYhszAYzTwC8sTFqc7jY9wuA0CIyDBgDvII2IXMwmR8vQQAAAAAgGVo5gEAAAAAsAzNPAJy8vOVc9/dyrnvbq7iCViIDAPmII+AXcgsTMY586gUJyfH7RIAVAEZBsxBHgG7kFmYiiPzAAAAAABYhmYeAAAAAADL0MwDAAAAAGAZmnkAAAAAACxDMw8AAAAAgGW4mj0Ci4pS1Bmt/MsALEOGAXOQR8AuZBYGo5lHQJ64OCX/8Rm3ywAQIjIMmIM8AnYhszAZXy8BAAAAAGAZjswDtZBn3MSQH+ukp4WxEgAAAACh4Mg8AnIKCpTzcJpyHk6TU1DgdjkAgkSGAXOQR8AuZBYm48g8AnMcOQcP+pdhhqocXUctQ4YBc5BHwC5kFgbjyDwAAAAAAJahmQcAAAAAwDI08wAAAAAAWIZmHgAAAAAAy9DMAwAAAABgGa5mj8A8HkU1a+ZfBmAZMgyYgzwCdiGzMBjNPALyxMcr+blxbpcBIERkGDAHeQTsQmZhMqbZAwAAAABgGZp5AAAAAAAsQzOPgJyCAh0bna5jo9PlFBS4XQ6AIJFhwBzkEbALmYXJOGcegTmOfNnZ/mUAliHDgDnII2AXMguDcWQeAAAAAADL0MwDAAAAAGAZmnkAAAAAACxDMw8AAAAAgGVo5gEAAAAAsAxXs0dgHo88p5ziXwZgGTIMmIM8AnYhszAYzTwC8sTHK+WliW6XASBEZBgwB3kE7EJmYTKm2QMAAAAAYBmaeQAAAAAALEMzj4CcwkIdG/O4jo15XE5hodvlAAgSGQbMQR4Bu5BZmIxz5hGYzyffrp3+ZdRunnGhnzfmpKeFsRJUGhkGzEEeAbuQWRiMI/MAAAAAAFiGZh4AAAAAAMvQzAMAAAAAYBmaeQAAAAAALEMzDwAAAACAZbiaPSrFk5LidgkAqoAMA+Ygj4BdyCxMRTOPgDwJCUp55TW3ywAQIjIMmIM8AnYhszAZ0+wBAAAAALAMzTwAAAAAAJZhmj0CcgoLlfv8c5KkpJGj5YmLc7kiAMEgw4A5yCNgFzILk9HMIzCfTyVbv/IvA7AMGQbMQR4Bu5BZGIxp9gAAAAAAWIZmHgAAAAAAy9DMAwAAAABgGZp5AAAAAAAsQzMPAAAAAIBluJo9Kicu3u0KAFQFGQbMQR4Bu5BZGIpmHgF5EhJUd9oMt8sAECIyDJiDPAJ2IbMwGdPsAQAAAACwDM08AAAAAACWYZo9AnIKC5U34SVJUuIDD8sTF+dyRQCCQYYBc5BHwC5kFiajmUdgPp+K/7XJvwzAMmQYMAd5BOxCZmEwptkDAAAAAGAZjswDiBjPuIkhP9ZJTwtjJQAAAIDdODIPAAAAAIBlaOYBAAAAALAMzTwAAAAAAJahmQcAAAAAwDJcAA8BeRISVPetd9wuo8aqykXhgMogw4A5yCNgFzILk3FkHgAAAAAAy9DMAwAAAABgGabZIyCnsFB5r74iSUocfp88cXEuVwQgGGQYMAd5BOxCZmEyjswjMJ9PxZnrVJy5TvL53K4GQLDIMGAO8gjYhczCYDTzAAAAAABYhmYeAAAAAADL0MwDAAAAAGAZmnkE5dnnn1fnzp2VkpKixo0b6+qrr1ZWVtavPubjjz+Wx+Mp92/r1q3+v/nwww/l9XpVr149DRkyRIWFhf77jhw5Iq/Xqz179lTbegG1SU5Ojh566CG1bNlSiYmJ6t69uzIzM/33V5RXj8ej559//qTPOWPGjAofk5+f7/+bv/71r2revLkaNGigkSNHlnn8N998I6/Xq6NHj4Z/hQFDsU0F7MM2FCahmUdQVq5apREjRujTTz/Vhx9+qOLiYvXp00c///xzwMdmZWXpu+++8/9r06aNJMnn82nQoEEaPny41q5dq4yMDE2dOtX/uEceeUTDhw9XixYtqm29gNpk2LBh+vDDD/XWW2/piy++UJ8+fXTZZZcpOztbksrk9LvvvtMbb7whj8ej66677left27duuUem5CQIEk6ePCghg0bpnHjxmnJkiWaOXOmFi1a5H/svffeq+eee05169atvhUHDMM2FbAP21CYhJ+mQ1A+WLhQnn9/sEjS9OnT1bhxY61fv149e/b81cc2btxY//Vf/1Xu9oMHD+rAgQO67777lJCQoIEDB2rLli2SpDVr1uizzz7T5MmTw7oe4eYZN9HtEoBKycvL07x587RgwQJ/Zp988kn9/e9/11/+8heNHTtWTZo0KfOYBQsWqFevXmrVqtWvPrfH4yn32FI7d+5UvXr1dNNNN0mSevXqpS1btqh///6aNWuW4uLidO2114ZhDQF7sE0F7MI2FKbhyDwCi49XytTpSpk6XYqPL3PXkSNHJEkNGjQI+DQdO3ZU06ZN1bt3by1fvtx/e6NGjdS0aVMtXbpUeXl5WrVqlTp06KDCwkLde++9evXVVxUdHR3edQJqkxMyXBwdrZKSEv+3/aUSExO1evXqcg/94YcftGjRIt15550BX+bYsWNq2bKlfvOb32jAgAHauHGj/742bdooNzdXGzdu1OHDh5WZmakOHTro8OHDGjNmjCZNmlT19QRswDYVsAvbUBiMZh4BeTweeRISjv/zePy3O46j3/3ud+rRo4fat29/0sc3bdpUr732mubNm6f58+crNTVVvXv31sqVK/3PP2fOHD399NNq166dOnbsqKFDh+q5555T7969lZiYqIsuukipqal8WAEhODHDdevW1YUXXqinn35a3377rUpKSvT2229r3bp1+u6778o9dubMmUpJSQn4jX/btm01Y8YMLVy4UO+8844SEhJ00UUX6euvv5Yk1a9fXzNnztTgwYPVpUsXDR48WH379lV6errS0tK0a9cudezYUe3bt9d7771XLeMAmIBtKmAXtqEwmcdxHCfQH+Xm5mrIkCGaOXOmkpKSIlEXLDBixAgtWrRIq1ev1m9+85ugHnvllVfK4/Fo4cKFFd6/bds29e/fXxs3blTPnj310EMP6YorrlD79u21bNkydejQIRyrEDZMs69+Tnqa2yXUGDt27NDQoUO1cuVKRUdHq1OnTvJ6vdqwYYN/Om6ptm3b6vLLL9fEicG9x30+nzp16qSePXtqwoQJFf7Nxx9/rJEjR2rFihVq3bq13nnnHTVp0kRdunTR119/rcaNG4e8jiiL7bj52KYCdmAbCpNwZB4BOUVFypvyF+VN+YucoiJJUlpamhYuXKjly5cHvdMhSd26dfN/21ju9RxHd999t1544QX5fD5t3LhR119/vRo3bqyLL75YK1asqNL6ALXNf2b4zDPP1IoVK3Ts2DHt3btXGRkZKioq0hlnnFHmcatWrVJWVpaGDRsW9GtGRUWpc+fOJ815QUGB7rvvPk2ZMkXbt29XcXGxLr74YqWmpsrr9WrdunUhrStgOrapgF3YhsJkNPMIrKRERatXqmj1SjnFxbr//vs1f/58/d///V+5D67K2rhxo5o2bVrhfdOmTVPDhg01cOBAlZSUSJKK/r3DU1RU5L8NQCWdkGGdkJ86deqoadOm+vHHH7VkyRJdddVVZR42bdo0nX/++Tr33HODfknHcbRp06aT5vzpp5/Wf//3f6tTp04qKSlRcXGx/z5yjhqNbSpgF7ahMBhXs0dQRjz0kN6ZPVsLFixQSkqKvv/+e0lSvXr1lJiYKEl69NFHlZ2drTfffFOSNH78eJ1++uk6++yzVVhYqLffflvz5s3TvHnzyj3//v37NXbsWK1Zs0bS8XOEzjrrLI0fP159+vTRRx99pMceeyxCawvUTEuWLJHjOEpNTdX27ds1cuRIpaam6re//a3/b44ePaq5c+fqhRdeqPA5Bg8erGbNmunZZ5+VJD311FPq1q2b2rRpo6NHj2rChAnatGlThVfN/vLLLzV79mxt2rRJ0vFpiFFRUZo2bZqaNGmirVu3qnPnzuFfccAwbFMB+7ANhUlo5hGUV197TZJ0ySWXlLl9+vTpuuOOOyQd/33NPXv2+O8rLCxUenq6srOzlZiYqLPPPluLFi1Sv379yj3/gw8+qPT0dDVr1sx/24wZMzRkyBBNmDBBI0eOVJcuXcK/YkAtcuTIET366KPat2+fGjRooOuuu07PPPOMYmNj/X/z7rvvynEc3XLLLRU+x549exQV9cvkrp9++kl33323vv/+e9WrV08dO3bUypUry+W1dMrvSy+9pDp16kg6fhXgGTNmaMSIESooKNCkSZPKfAYANRXbVMA+bENhEi6Ah4Cc/Hzl3HX828aUqdPL/CYujuMCeNWPC+CFjgyD7bg5yCNgFzILk3HOPAAAAAAAlqGZBwAAAADAMjTzAAAAAABYhgvgIbD4eCVPnuJfBmAZMgyYgzwCdiGzMBjNPALyeDzy1K3rdhkAQkSGAXOQR8AuZBYmY5o9AAAAAACW4cg8AnKKipT/17ckSQmDbpfnhN/RBGA+MgyYgzwCdiGzMBlH5hFYSYmKPvpQRR99KJWUuF0NgGCRYcAc5BGwC5mFwWjmAQAAAACwDM08AAAAAACWoZkHAAAAAMAyNPMAAAAAAFiGZh4AAAAAAMvQzAMAAAAAYBl+Zx6BxcUp+cUJ/mUAliHDgDnII2AXMguD0cwjIE9UlDyNGrldBoAQkWHAHOQRsAuZhcmYZg8AAAAAgGU4Mo+AnOJiFcydLUmKv+EmeWJ42wA2IcOAOcgjYBcyC5NxZB6BFRercPH7Klz8vlRc7HY1AIJFhgFzkEfALmQWBqOZBwAAAADAMjTzAAAAAABYhmYeAAAAAADL0MwDAAAAAGAZmnkAAAAAACxDMw8AAAAAgGX4oUQEFhenOs/+2b8MuMEzbmJIj3PS08JciYXIMGAO8gjYhczCYDTzCMgTFaXo3zR3uwwAISLDgDnII2AXMguTMc0eAAAAAADLcGQeATnFxSpY+HdJUvzAq+WJ4W0D2IQMA+Ygj4BdyCxMxrsRgRUXq/Bv8yRJ8f0GSHyIAXYhw4A5yCNgFzILgzHNHgAAAAAAy/DVEvBvoV4tHQAAAAAijSPzAAAAAABYhmYeAAAAAADL0MwDAAAAAGAZmnkAAAAAACzDBfAQWFyc6jw11r8MwDJkGDAHeQTsQmZhMJp5BOSJilJ0qzPdLgNAiMgwYA7yCNiFzMJkTLMHAAAAAMAyHJlHQE5xsQqXfCBJiuv73/LE8LYBbEKGAXOQR8AuZBYm492IwIqLVfDuLElSXO/LJT7EALuQYcAc5BGwC5mFwZhmDwAAAACAZWjmAQAAAACwDM08AAAAAACWoZkHAAAAAMAyNPMAAAAAAFiGZh4AAAAAAMvw2woILC5OSY/9wb8MwDJkGDAHeQTsQmZhMJp5BOSJilLMWe3cLgNAiMgwYA7yCNiFzMJkTLMHAAAAAMAyHJlHQE5xsYqW/58kKbbXpfLE8LYBbEKGAXOQR8AuZBYm492IwIqLlf/mdElS7P/rKfEhBtiFDAPmII+AXcgsDMY0ewAAAAAALEMzDwAAAACAZWjmAQAAAACwDM08AAAAAACWoZkHAAAAAMAyNPMAAAAAAFiG31ZAYLGxSvyfkf5lAJYhw4A5yCNgFzILg9HMIyBPdLRiz+vkdhkAQkSGAXOQR8AuZBYmY5o9AAAAAACW4cg8AnKKi1W0do0kKbb7RfLE8LYBbEKGAXOQR8AuZBYm492IwIqLlT/1VUlSbJeuEh9igF3IMGAO8gjYhczCYEyzBwAAAADAMjTzAAAAAABYhmYeAAAAAADL0MwDAAAAAGAZmnkAAAAAACxDMw8AAAAAgGX4bQUEFhurxPsf9C8DsAwZBsxBHgG7kFkYjGYeAXmioxXbtZvbZQAIERkGzEEeAbuQWZiMafYAAAAAAFiGI/MIyCkpUfFnmZKkmAs6yxMd7XJFAIJBhgFzkEfALmQWJqOZR2BFRcqb9LIkKWXqdIkPMcAuZBgwB3kE7EJmYTCm2QMAAAAAYBmaeQAAAAAALEMzDwAAAACAZWjmAQAAAACwDBfAA1CjecZNDPmxTnpaGCsBAAAAwocj8wAAAAAAWIYj8wgsJkYJdw33LwOwDBkGzEEeAbuQWRiMdyQC8sTEKK7nxW6XASBEZBgwB3kE7EJmYTKm2QMAAAAAYBmOzCMgp6RExV/8S5IUc8658kRHu1wRgGCQYcAc5BGwC5mFyWjmEVhRkfJeeF6SlDJ1usSHGGAXMgyYgzwCdiGzMBjT7AEAAAAAsAzNPAAAAAAAlmGaPWoUz7iJbpcAAAAAANWOI/MAAAAAAFiGZh4AAAAAAMvQzAMAAAAAYBnOmUdgMTFKGPxb/zIAy5BhwBzkEbALmYXBeEciIE9MjOIu7+N2GQBCRIYBc5BHwC5kFiZjmj0AAAAAAJbhyDwCcnw+lWRtlSRFp7aVJ4rvgACbkGHAHOQRsAuZhclo5hFYYaFy//dpSVLK1OlSQoLLBQEIChkGzEEeAbuQWRiMr5YAAAAAALAMzTwAAAAAAJahmQcAAAAAwDI08wAAAAAAWIZmHgAAAAAAy9DMAwAAAABgGX6aDoHFxCj+5lv9ywAsQ4YBc5BHwC5kFgbjHYmAPDExiu9/pdtlAAgRGQbMQR4Bu5BZmIxp9gAAAAAAWIYj8wjI8fnk+2aXJCnq9DPkieI7IMAmZBgwB3kE7EJmYTKaeQRWWKifn/i9JCll6nQpIcHlggAEhQwD5iCPgF3ILAzGV0sAAAAAAFiGZh4AAAAAAMvQzAMAAAAAYBmaeQAAAAAALEMzDwAAAACAZbiaPQCchGfcxJAf66SnhbESAAAAoCyaeQQWE6O4a67zLwOwDBkGzEEeAbuQWRiMdyQC8sTEKOHa690uA0CIyDBgDvII2IXMwmScMw8AAAAAgGU4Mo+AHJ9Pvm+zJUlRpzWTJ4rvgACbkGHAHOQRsAuZhclo5hFYYaF+fnSUJCll6nQpIcHlggAEhQwD5iCPgF3ILAzGV0sAAAAAAFiGZh4AAAAAAMvQzAMAAAAAYBmaeQAAAAAALEMzDwAAAACAZWjmAQAAAACwDD9Nh8BiYhTXb4B/GYBlyDBgDvII2IXMwmC8IxGQJyZGCbcMcrsMACEiw4A5yCNgFzILkzHNHgAAAAAAy3BkHgE5Pp+cQ4ckSZ6GDeWJ4jsgwCZkGDAHeQTsQmZhMpp5BFZYqGO/e0CSlDJ1upSQ4HJBAIJChgFzkEfALmQWBqOZh3E84ya6XQIAAAAAGI15IgAAAAAAWIZmHgAAAAAAy9DMAwAAAABgGZp5AAAAAAAsQzMPAAAAAIBluJo9AouOVmzvy/3LACxDhgFzkEfALmQWBqOZR0Ce2Fgl3jHU7TIAhIgMA+Ygj4BdyCxMxjR7AAAAAAAsw5F5BOQ4jpycHEmSJyVFHo/H5YoABIMMA+Ygj4BdyCxMRjOPwAoKdGzEPZKklKnTpYQElwsCEBQyDJiDPAJ2IbMwGNPsAQAAAACwDM08AAAAAACWoZkHAAAAAMAyNPMAAAAAAFiGZh4AAAAAAMtwNXsAqAaecRNDfqyTnhbGSgAAAFAT0cwjsOhoxfbo6V8GYBkyDJiDPAJ2IbMwGM08AvLExirxnnvdLgNAiMgwYA7yCNiFzMJknDMPAAAAAIBlODKPgBzHkQoKjv9HfLw8Ho+7BQEIChkGzEEeAbuQWZiMI/MIrKBAOXf9Vjl3/faXDzMA9iDDgDnII2AXMguD0cwDAAAAAGAZmnkAAAAAACxDMw8AAAAAgGVo5gEAAAAAsAzNPAAAAAAAlqGZBwAAAADAMvzOPAKLilJM567+ZQCWIcOAOcgjYBcyC4PRzCMgT1yckh54yO0yAISIDAPmII+AXcgsTMbXSwAAAAAAWIZmHgAAAAAAy9DMIyAnP19Hb79FR2+/RU5+vtvlAAgSGQbMQR4Bu5BZmIxmHgAAAAAAy9DMAwAAAABgGa5mj2rjGTfR7RIAAAAAoEbiyDwAAAAAAJahmQcAAAAAwDI08wAAAAAAWIZz5hFYVJRizj3PvwzAMmQYMAd5BOxCZmEwmnkE5ImLU1L6I26XASBEZBgwB3kE7EJmYTK+XgIAAAAAwDI08wAAAAAAWIZmHgE5+fk6eucdOnrnHXLy890uB0CQyDBgDvII2IXMwmScM4/KKSxwuwIAVUGGAXOQR8AuZBaG4sg8AAAAAACWoZkHAAAAAMAyNPMAAAAAAFiGZh4AAAAAAMtwATwAMIxn3MSQH+ukp4WxEgAAAJiKZh6BRUUpuu1Z/mUAliHDgDnII2AXMguD0cwjIE9cnOo8PsbtMgCEiAwD5iCPgF3ILEzG10sAAAAAAFiGZh4AAAAAAMswzR4BOfn5Ova7ByRJyS9OkCchweWKAASDDAPmII9m4sKjOBkyC5PRzKNSnJwct0sAUAVkGDAHeaxZ+CKg5iOzMBXT7AEAAAAAsAzNPAAAAAAAlqGZBwAAAADAMpwzDwAAAOtV5dx1ALARR+YBAAAAALAMR+YRWFSUos5o5V8GYBkyDJiDPAJ2IbMwGM08AvLExSn5j8+4XQaAEJFhwBzkMTCmy8MkZBYmo5nHr2KDCgAAAADmYa4IAAAAAACWoZlHQIm+En2+Zb0+37Jeib4St8sBECSnoEA5D6cp5+E0OQUFbpcD1GrkEbALmYXJmGaPgDyO1LKowL8MwDKOI+fgQf8yABeRR5ygKqczOulpYawEJ0VmYTCOzAMAAAAAYBmOzAMAAACWCfWoPkf0gZqDI/MAAAAAAFiGZh4AAAAAAMvQzAMAAAAAYBnOmUdAjkf6Kj7RvwzAMh6Popo18y8DcBF5BOxCZmEwmnkElBcVrW5tO7pdBoAQeeLjlfzcOLfLACDyCNiGzMJkNPMAAAAIq6r8fjoAoHI4Zx4AAAAAAMtwZB4BJfpKtHzb55KkXt4OyouKdrkiAMFwCgr08xOPS5LqPPWMPPHxLlcE1F7kEW6ryqyJ2vgb9WQWJqOZR0AeRzqrIM+/DMAyjiNfdrZ/GYCLyCNgFzILgzHNHgAAAAAAy9DMAwAAAABgGZp5AAAAAAAswznzAFCDVHRho6SSEn337+Xkl19VbnT5i1jWxosaAQAA2IxmHgAAAOXwW/EAYDaaeQTkeKTdsfH+ZQB2IcOAQTweeU45xb8MwHBkFgajmUdAeVHR6tDufLfLABAiMgyYwxMfr5SXOOIN2ILMwmRcAA8AAAAAAMtwZL4W4Jw3AAAAVFVV9im50CoQfhyZR0AJvhIt3/YvLd/2LyX4StwuB0CQyDBgDqewUMfGPK5jYx6XU1jodjkAAiCzMBlH5hFQlCN1yvvZvwzALmQYMIjPJ9+unf5lAIYjszAYzTwAgKmTAAAAlmGaPQAAAAAAlqGZBwAAAADAMkyzBwAAqKH4RRsAqLk4Mg8AAAAAgGU4Mo9KORjNWwWwGRkGzOFJSXG7BABBILMwFXt3CCg3Olpntu/idhkAQkSGAXN4EhKU8sprbpcBoJLILExGMw8AAACgWvETqED4cc48AAAAAACW4cg8Akrwlei9nV9Jkq5vdZbyo6JdrghAMMgwYA6nsFC5zz8nSUoaOVqeuDiXKwLwa8gsTEYzj4CiHOn//XzUvwzALtWdYaZOAkHw+VSy9Sv/MgDDkVkYjGbeEvxOLAAAtRP7AACAitDMAwAAADAWM8CAinEBPAAAAAAALEMzDwAAAACAZWjmAQAAAACwDOfMo1J+juJ7H8BmZBgwSFy82xUACAaZhaFo5hFQbnS0Tjunm9tlAAgRGQbM4UlIUN1pM9wuA6g1Qr14XumF88gsTEYzDwAAEAH8xBwAIJxo5iOIjTgAlMXPDQEAAISGZh4Bxft8euubrZKk209vqwLOvQWsQoYBc5BHwC5OYaHyJrwkSUp84GF54uJcrgj4Bc08Aop2HPXN+cm/DMAuZBgIn6rOsiOPgGV8PhX/a5N/GTAJzTwAwEpM0QcAVJfSbUxSSYm++/dtyS+/qtzo6ICPZRuDSKGZDwHnvgMAAAAA3FRrm3kacgAAaif2AQBUJ2aOIVK46goAAAAAAJaptUfmAQAAAMAkHNVHMCrVzDv/vtpqXl5e2AuoN2FK2J+zMvgWo/JifCXK+/d7IKakWDFcfRewChkuL/a5lyL+mkceuKfanjsxMVEej+ek91fndtwtVdl/cHMfgDwCdrEpszVt21bbBNqWV8TjOIHfkYcOHdLw4cNDLgwAAFSfmTNnKikp6aT3sx0HAMBsgbblFalUM+/z+fTjjz8qISEh6G8LaoOff/5ZPXv21MqVK1WnTh23y7ECYxYcxit4jFnwGLPgmDRegb7Nd3M7btI4mY6xqhzGqfIYq8phnCqPsaqcUMYplCPzlZppFhUVpYYNGwb1xLWJz+eTz+dTYmJi0N+m1FaMWXAYr+AxZsFjzIJj03i5uR23aZzcxlhVDuNUeYxV5TBOlcdYVU6kxomr2QMAAAAAYBmaeQAAAAAALEMzHwZxcXG6//77FRcX53Yp1mDMgsN4BY8xCx5jFhzGq3IYp8pjrCqHcao8xqpyGKfKY6wqJ1LjVKkL4AEAAAAAAHNwZB4AAAAAAMvQzAMAAAAAYBmaeQAAAAAALEMzDwAAAACAZWjmAQAAAACwTIzbBbjt888/18SJE7Vp0yYVFRWpdevWGjJkiK688sqQnq+oqEjXX3+9tm7dqjPOOEP//Oc/I/K6kRLp8frhhx/0wQcfaOXKldq5c6cOHjyoevXqqVOnTho2bJjOPffccKxWtXLrPXaiqVOnaty4cZKk2bNn67zzzgvptSPFzTH78MMPNWvWLG3ZskV5eXk65ZRTdN5552nkyJFq2rRpqKtUrdwYL8dx9OGHH+qtt97Srl27lJOToyZNmqhr166666671Lx586quVrWq6pitW7dOgwcPPun9J8uZrZ/9J1qwYIHWr1+vzZs3a9u2bSoqKtKzzz6ra6+9Nqjn8fl8mjVrlmbPnq3du3crKSlJXbt21cMPP6zTTz+9eoqPoHCM06FDh/Tee+/pyy+/1ObNm5WdnS1JysrKqq6yXRGOsfrss8+0bNkyZWRkKDs7W7m5uWrWrJl69+6te+65R3Xr1q3GNYiMcIzTunXrNGfOHG3ZskUHDhxQUVGRmjRpok6dOumuu+5Sq1atqnENIidcn1MnCmX/y3Thek+Fsj20TTjfU8eOHdMbb7yhpUuXau/evYqNjVXz5s3Vu3dv3X///UE9V61u5tetW6c777xTsbGx6t+/v1JSUrR06VKlp6crOztbw4cPD/o5X3nlFe3ZsyfirxsJbozXW2+9palTp6pFixbq3r27GjZsqN27d2vZsmVatmyZXnjhBfXr168qq1Wt3HqPnWjHjh2aMGGCkpKSlJubG/TrRZpbY+Y4jp544gnNnj1bLVq0UL9+/VSnTh3t379fmZmZys7ONrKZd2u8/vSnP2n69Olq1KiRevfureTkZG3dulVz5szR+++/r3fffVderzfU1apW4RyzLl26qEuXLuVub9KkSbW+rptefvllZWdnq379+mrcuLG/wQzWE088oTlz5qh169a67bbbdOjQIS1evFhr1qzRu+++q9atW4e58sgKxzht375dL774ojwej1q2bKnExETl5eVVQ7XuCsdYPfjgg/rxxx91/vnn66qrrpLH41FGRoZef/11LV26VO+++64aNmxYDdVHTjjGae3atVq/fr06dOigHj16KDY2Vjt37tSCBQv0/vvva+rUqerWrVs1VB9Z4fqcOlGw+182COc4BbM9tFG4xurbb7/VkCFDtHfvXnXv3l0XX3yxCgsLtWfPHi1ZsiToZl5OLVVUVORcdtllTvv27Z0vv/zSf3tOTo7Tv39/p127ds6uXbuCes7Nmzc77dq1c958803H6/U6ffv2jcjrRoJb47VkyRInMzOz3O2ZmZnO2Wef7XTp0sUpKCgIen0iwa0xO1FxcbFz3XXXOddff72Tnp7ueL1eZ+PGjSGsTWS4OWYzZ850vF6v89RTTznFxcUV1mYat8Zr//79Ttu2bZ1evXo5OTk5Ze6bPn264/V6ndGjR4e0TtUtXGP26aefOl6v15kwYUJEX9cEa9ascfbt2+c4juNMmTLF8Xq9zrx584J6jk8++cTxer3OrbfeWuYzfO3atU5qaqozaNCgsNbshnCM04EDB5yMjAx/zvr27et4vd6w1+q2cIzVlClTnB9++KHMbT6fz3niiSccr9frPPnkk2Gr1y3hGKf8/PwKb1+7dq3j9Xqda6+9tsp1miAcY3WiYPe/bBGOcQp2e2ircIxV6X55hw4dnE8++aTc/aHsa9bac+Y//fRT7dmzRwMGDFC7du38tycnJ+u+++5TcXGx5s+fX+nnKyws1OjRo3Xuuefqtttui9jrRopb49WnTx9dcMEF5W6/4IIL1LVrV/3000/GTjd0a8xONHXqVG3dulX/+7//q+jo6KDXIdLcGrP8/HxNnjxZzZs312OPPVbhWMXEmDeRya3xys7Ols/nU6dOnZScnFzmvksuuUSSdPjw4eBWJkLc+gy29bO/It27d1ezZs2q9Bxz586VJD300EOKi4vz337hhReqR48eyszM1K5du6r0Gm4Lxzidcsop6ty5c7mc1TThGKu7775bjRs3LnObx+PRfffdJ0nKzMys0vObIBzjFB8fX+HtF154oerVq1djjjyHY6xKhbL/ZYtwjlNNF46xWrJkib744gsNHTq0whkwoexrmrd3GiEZGRmSpB49epS776KLLirzN5UxadIk7d69WwsWLJDH44nY60aKW+P1a0rf8CY2WZL7Y7Zt2zZNmjRJ9957r9q0aVPp13GTW2O2Zs0a/fTTT7rmmmvk8/m0dOlSffPNN0pJSVH37t3VsmXLINckMtwar5YtWyo2NlYbNmzQsWPHyjQaK1askCRjp2mGe8y++eYbvfnmm8rPz9dpp52m7t27q0GDBtX+urZbt26dkpKS1KlTp3L39ejRQ6tWrVJmZqbOOOMMF6pDTVK6j2DDF9pu2rhxo44cOaLzzz/f7VKME4591tqgstvD2mzx4sWSpCuuuELfffedPv74Y+Xk5Kh58+bq2bOn6tSpE/RzmtkFRcA333wjSRXupNerV0/169fX7t27K/Vcn3/+uV5//XU9/PDDAXc8wvm6keTWeJ3Mt99+q7Vr16pRo0bGnpfr5pgVFxdr9OjROvPMM3X33XcHVbeb3BqzzZs3Szq+szdw4MAyRwSjoqJ0xx136JFHHqnkWkSOW+NVv359Pfzww/rzn/+sfv366dJLL1WdOnW0bds2ffLJJ7rpppuMPXoR7s/g999/X++//77/vxMSEpSWlqZhw4ZV6+vaLDc3VwcOHJDX662wwSq9+F3pmAFVMW/ePEm/fGmG49atW6eMjAwVFhZq9+7dWr58uerXr69HH33U7dKMEo591tqistvD2qx0f3P9+vV69tlnVVhY6L+vQYMGGj9+vLp27RrUc9baZv7YsWOSpJSUlArvT05O1vfffx/weQoLC/Xoo4/qrLPO0tChQyP2upHm1nhVpKioSKNGjVJhYaHS09ON/bbdzTF79dVXlZWVpTlz5ig2NrbyRbvMrTE7dOiQJGn69Olq166d5s6dqzPPPFNfffWV/vCHP+iNN95Q8+bNdeuttwaxNtXPzffYnXfeqcaNG2vMmDF65513/Ld37NhRAwcONPZ9F64xa9CggUaNGqVLLrlEp512mo4ePap169Zp3Lhxev7555WcnKybb7457K9bE+Tk5EjSSaeOl95eOmZAqL766itNnjxZDRs2pKH4DxkZGZo0aZL/v1u2bKkXX3xR7du3d7Eqs4Rjn7U2CHZ7WJuV7m+OHTtWQ4cO1W233aa4uDgtWrRIf/rTnzRixAgtXry43ClDv6bWnjMfLuPHj9fu3butOSfZbVUdL5/Pp8cee0yZmZm68cYbdfXVV4e/SMMEO2Zbt27Vq6++qqFDh+rss8+OQIXmCXbMHMeRJMXGxmry5Mnq0KGD6tSpowsuuEATJkxQVFSUpk+fXt1luyaUXL7yyit69NFHdc8992jFihXauHGjZs2apZKSEg0ePFhLly6t5qrd1aZNG915550688wzlZiYqFNPPVUDBw7U66+/rtjYWE2cOFE+n8/tMoFaa+/evbrnnntUUlKiF198kem+/yEtLU1ZWVnauHGj5s6dq1atWumWW27RP/7xD7dLMwb7+JXD9rDySvc3L7nkEqWnp6tJkyZq0KCBbr/9dt1xxx3KycnRe++9F9Rz1tpmvvSb/9IjBP/p2LFjJz2CUurLL7/UjBkzNHz4cKWmpkbsdd3g1nidyHEc/f73v9fChQs1cOBAPfXUU0E/RyS5NWaPPPKImjdvrrS0tOAKNoDbuWzfvr1OPfXUMve1adNGzZs31549e3T06NFKPV+kuDVen3zyiV5++WUNGjRIw4cPV5MmTZSUlKTzzz9fU6ZMUXx8vJ599tngViZCqvsz2Ov16txzz9XBgwfLTJu39bO/OpSu58mOvJfeXtMv+obqk52drSFDhujw4cOaMGGCsdfwMEFSUpI6dOigSZMmqVWrVhozZoyxFzCNpKrus+Lk28ParHS7dumll5a7r1evXpJ+mYpfWbW2mS89J6+iN9eRI0f0448/BrzoVVZWlkpKSjRx4kSlpqaW+SdJu3btUmpqapmrsYfjdd3g1niVKj0iP2/ePA0YMEDPPfecoqLMfvu6NWZbt27Vzp07dc4555T5+7/97W+SpJtuukmpqalatmxZmNY0fNwas1atWkk6+RTo0tvz8/ODXqfq5NZ4lV7krqLzuho0aKDU1FR9++23Ru4QRuIzuH79+pLKvl9s/eyvDklJSWrUqJH27dunkpKScveXnitfOmZAMPbt26fbb79d+/fv1/jx4/07yPh1MTEx6tq1q3Jzc/XFF1+4XY7rqrLPil9UtD2szUqvu1C3bt1y95XeVlBQENRz1tpz5jt37qwpU6Zo9erV6t+/f5n71qxZI0nq0qXLrz7H6aefruuvv77C+9577z2lpKSob9++SkxMDOvrusGt8ZKON/KPP/645s+fr379+unPf/6zFdOd3Bqzk/39Z599pm+++UaXXnqpGjRoYORPkbg1ZqVN6c6dO8s9pqioSHv27FFSUpJx0zTdGq+ioiJJJ//5udLbT/zJMVNU92dwcXGxtmzZIo/Ho6ZNm0bsdW3TpUsXLVq0SBs2bFDnzp3L3Ld69WpJKnc7EMi+ffs0ePBg7d+/Xy+99JIuu+wyt0uyyv79+yWZ+ytBkRTqPit+cbLtYW3WrVs3bdiwQdu3b1efPn3K3Ld9+3ZJCn7/POhfpq8hioqKnN69ezvt27d3tmzZ4r89JyfH6d+/v9OuXTtn586d/tsPHTrkbN++3Tl06FClnt/r9Tp9+/at8uuawq3xKikpcUaPHu14vV7ngQcecIqKiqq+MhHi1pidzCOPPOJ4vV5n48aNlX5MpLk5ZkOHDnW8Xq8zZ86cMrdPmjTJ8Xq9Tnp6eghrVL3cGq/333/f8Xq9Tv/+/Z2jR4+WuW/+/PmO1+t1rrnmmhDXqnqFa8w2bNjg+Hy+cs/9zDPPOF6v17nzzjur9Lq2mDJliuP1ep158+ZVeP/Jxu+TTz5xvF6vc+uttzoFBQX+29euXeukpqY6gwYNqta6Iy3UcfpPffv2dbxeb3WUaIxQx2rv3r1Or169nHbt2jlLliyJRKmuCnWcMjIyyn12OY7jrFq1yjn77LOd888/3/n555+rpWa3hCt/pYLd/7JFqOMU7PawJgh1rPbs2eO0b9/eufDCC53vv//ef3tOTo5z1VVXOV6v11m7dm1QtdTar95iYmI0duxYDRs2TLfeeqsGDBig5ORkLV26VPv27dNDDz1U5ico/vrXv2rSpEm6//77q3QucrCvawq3xmvy5MmaP3++kpKSdPrpp+svf/lLub+57LLLdNZZZ4X8GtXFrTGzmZtj9sQTT+jmm2/W73//ey1btkytWrXSli1b9Omnn6pZs2YaNWpUVVcv7NwaryuuuELvvvuuMjIy1KdPH1166aWqW7eusrKytGbNGsXFxemxxx4LxyqGXbjG7H/+538kHb96/6mnnqqcnBxlZmZq165dOu2008pd08PWz/6KzJ07V+vXr5ckbdu2zX9bRkaGpOOfyaVHRE82ft26ddMNN9yguXPn6pprrtHFF1+sQ4cOafHixUpOTtaTTz4Z2ZWqBuEYJ0kaPXq0f/nAgQPlbhs1apRxs4aCFY6xGjx4sLKzs3XeeecpKytLWVlZ5V7H9m1rOMbp3nvvVf369XXOOeeoSZMmKigoUFZWljIzMxUbG6uxY8cqKSkpwmsWfuHKX00XjnEKdntoq3CMVfPmzTVq1CiNHTtWAwcO1OWXX664uDh9/PHHys7O1k033aQLL7wwqLpqbTMvHd+ZmDVrliZMmKAPPvhARUVFat26tR588EENHDiwxr1uVblRd3Z2tqTjv0v86quvVvg3zZo1M7KZl+z9f+0mt8asRYsWmjdvniZMmKBVq1ZpzZo1OuWUUzRo0CCNGDFCDRs2rLbXrgo3xis6OlrTpk3TzJkz9cEHH2jRokUqKipSw4YNNWDAAN1zzz3yer3V8trhEI4xu/nmm7Vq1SplZGToxx9/VExMjFq0aKHhw4dr6NChqlevXrW8rgnWr1/vvwZHqQ0bNmjDhg2Sjn8mV2Z68x//+EelpqZq9uzZeuutt5SUlKRevXrVmN9zDtc4/edz/Odt999/v/XNfDjGqnR/YdOmTdq0aVOFf2N7oxaOcUpLS9OqVau0fv16HT582D8F+oYbbtCQIUPUpk2baqs/ksKVv5ouHOMUyvbQRuF6T91+++1q1qyZpk2bpkWLFqmkpEStW7fW8OHDdeONNwZdl8dx/n2NfAAAAAAAYAWzLwcOAAAAAADKoZkHAAAAAMAyNPMAAAAAAFiGZh4AAAAAAMvQzAMAAAAAYBmaeQAAAAAALEMzDwAAAACAZWjmAQAAAACwDM08AAAAAACWoZkHAAAAAMAyNPMAAAAAAFiGZh4AAAAAAMv8f0gUrVKLtxqHAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_, axes = plt.subplots(1, 2, figsize=(10, 4), sharey=True)\n",
    "for sample, title, ax in zip((theta, phi),\n",
    "                             (r\"$\\theta$\", r\"$\\phi$\"),\n",
    "                             axes.ravel()):\n",
    "    # Histogram plots with 30 bins for theta and phi\n",
    "    ax.hist(sample, bins=30, density=True)\n",
    "    # Compute 2.5% and 97.5% quantile approximation using samples\n",
    "    q_s = np.quantile(sample, [0.025, 0.975])\n",
    "    ax.axvline(q_s[0], color=\"C2\", ls=\"--\")\n",
    "    ax.axvline(q_s[1], color=\"C2\", ls=\"--\")\n",
    "    ax.text(q_s[0], 20, \"2.5%\", color=\"k\")\n",
    "    ax.text(q_s[1], 20, \"97.5%\", color=\"k\")\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Authors:\n",
    "- Aki Vehtari <aki.vehtari@aalto.fi>\n",
    "- Tuomas Sivula <tuomas.sivula@aalto.fi>\n",
    "- Osvaldo A. Martin <osvaldo.martin@aalto.fi>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bda_py_demos",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
